% Test passed on 08-10-12
% 

n = 2; t = (0:n)/n;
[V, T] = init_mesh_rectangle(t, t, 'type3C');
[T, E, ET, TE] = build_fem_mesh(V, T);

FE_Order = 3;

lambda = 100;
mu = 100;
R = 0.5;

[dof_map, n_dof] = distribute_dof(T, TE, ET, FE_Order, 0);  

[dof_x, dof_y] = build_dof_coordinate(dof_map, V, T, FE_Order);

bdr_Dirichlet = find(ET(:,2)==0);
[dof_boundary, ub_1, ub_2] = dirichlet_bb(dof_map, V, T, TE, ET, FE_Order, bdr_Dirichlet, @fun_u, lambda, mu, R);

[H, constrains] = constrains_C0(dof_map, TE, ET, FE_Order); 

% ZERO = zeros(size(V,1),1); 
% trisurf(T, V(:,1), V(:,2), ZERO);
% axis tight;view(0, 90);  
% hold on;
% text(dof_x(constrains(:,1)), dof_y(constrains(:,1))+ 0.05, num2str(constrains(:,1)));
% text(dof_x(constrains(:,2)), dof_y(constrains(:,2)), num2str(constrains(:,2)));
% hold off;

% plot the dof element-by-element
dofs_per_elem = (FE_Order + 1)*(FE_Order + 2)/2;
n_elem = size(T,1);
[I, J, K] = indices(FE_Order);
I = I/FE_Order;
J = J/FE_Order;
K= K/FE_Order;

ZERO = zeros(size(V,1),1); 
trisurf(T, V(:,1), V(:,2), ZERO);
axis([0,1,0,1]);axis tight; view(0, 90);  
hold on;

for k = 1:n_elem
    
    V1 = V(T(k,1),:);
    V2 = V(T(k,2),:);
    V3 = V(T(k,3),:);
    
    VC = (V1 + V2 + V3)/3;
    
    shrink = 0.8;
    V1 = VC + shrink*(V1 - VC);
    V2 = VC + shrink*(V2 - VC);
    V3 = VC + shrink*(V3 - VC);
    
    coord = I*V1 + J*V2 + K*V3;
    
%     dofs = (1:dofs_per_elem) + (k-1)*dofs_per_elem;
    dofs_str = num2str(dof_map(:,k));
    
    text(coord(:,1), coord(:,2), dofs_str);
end

hold off;


% then you can validate the constrains according to the figure